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ABSTRACT 

Wavenumber transition and hysteresis in a highly unstable baroclinic flow are investigated using a high- 
resolution spectral numerical model. As the flow becomes more supercritical, the dominant wave gradually shifts 
from the most unstable wave predicted by the linear theory to a longer wave with a larger time-averaged 
amplitude, while the rectified mean flow attains a stronger shear at the center of the channel. The numerical 
results display a complex hysteresis behavior, which occurs not only between the states of different dominant 
wavenumbers, but also between the states of identical dominant wavenumber but of different dynamic charac- 
teristics. In a certain parameter range three stable states, each with different dominant wavenumber, are possible, ) 

and in another parameter range four stable states are possible, among them three stable states with an identical j. 

dominant wave. The numerical results suggest that a multiple weather regime exists even without external forcing i 

in which the flow aperiodically varies between two distinct behaviors. The effects of stable higher harmonics J 

are assessed and it is found that their presence contributes not only to the better approximation of the model ) 

solutions but also to the selection of the final equilibrium state, due to the chaotic nature of the initial transient 
period. - 


1. Introduction 

Most theories concerning the finite-amplitude behav- 
ior of unstable baroclinic waves have usually focused 
on the dynamics of a single wave. The restriction to a 
single wave is justifiable in weakly nonlinear theories 
since, with the assumption of small supercriticality, the 
dynamics for a single unstable wave become closed in 
a consistent way if the initial wave spectrum is also 
limited to that wave (e.g., Pedlosky 1970; Drazin 1972; 
Chou and Loesch 1986a,b). However, these theories 
can only be applied to a narrow range in parameter 
space where the deviation from marginal instability is 
small. The single wave restriction is sometimes justi- 
fied since single wave states are frequently observed in 
annulus experiments, although it is unclear whether the 
observed wave represents the linearly most unstable 
wave. This wave scale probably dominates the wave 
spectrum at the initial growing stage, but there is no 
reason to believe that, as the zonal flow is altered by 
the growing unstable waves, it will remain dominant at 
finite amplitude. As the flow becomes more supercrit- 
ical, linear theories predict that the fastest growing 
wave may gradually shift to a shorter wave. In contrast, 
laboratory experiments and weakly nonlinear theories 
have shown that the long wave may grow at the ex- 
pense of the linearly most unstable wave and eventually 
become dominant at finite amplitude. 
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Another fundamental problem in geophysical-re- 
lated fluid dynamics is the presence of hysteresis. In 
general, hysteresis occurs when, for some points in pa- 
rameter space, different fluid behavior can be observed 
depending on the initial conditions of the experiment. 
Once the flow settles in a particular behavior, it remains 
in that state beyond the point in parameter space where 
it first developed. Hysteresis has been observed in the 
axisymmetric to wave regime transition (Fein 1973; 
Miller and Butler 1991). It has also been observed in- 
side the wave regime and is usually associated with 
wavenumber selection, where the wavenumber ob- 
served varies with initial condition. 

The presence of hysteresis and wavenumber transi- 
tion in a thermally driven rotating annulus has been 
shown by Buzyna et al. ( 1978), when the imposed tem- 
perature contrast is varied gradually. This study in- 
spired Pedlosky ( 1981 ) to analytically investigate the 
weakly nonlinear dynamics of wave ensembles in the 
limit of small Froude number, weak viscosity, and long 
waves, using a two-layer model on an /plane. His re- 
sults show that in the case of free unstable waves the 
wave realized at finite amplitude is not the linearly most 
unstable wave. Rather, a long wave, capable of achiev- 
ing the single largest steady amplitude, is favored in 
the competition for the potential energy of the basic 
state. The presence of conjugate waves, which are ca- 
pable of achieving identical final amplitude, implies the 
presence of hysteresis. Hart (1981) in his laboratory 
experiment containing two layers of immiscible fluid 
in a rotating frame shows that there is a succession of 
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Fig. I . Linear growth rate as a function of wavenumber k 
and vertical shear X for S = 0.1, r = 1, and r, =0.1. 


wavenumber transition to longer waves as the system 
becomes more supercritical. He constructed a simple 
analytical model that consists of only two adjacent un- 
stable waves closed to the intersection of their marginal 
curves to interpret his laboratory results. Through 
weakly nonlinear interactions between these waves and 
the mean flow, either one of these waves may become 
dominant at the expense of the other, or the waves may 
coexist. These results explained reasonably well the 
wavenumber transition in the immediate vicinity of the 
intersection point of two marginal curves for adjacent 
wavenumbers, but observations of annulus experiments 
showed that the wavenumber transition may not be lim- 
ited to that region. 

Recently, a numerical study of wavenumber selec- 
tion in a rotating, differentially heated annulus flow 
system was conducted by Lu et al. (1994) using a 
three-dimensional fully nonlinear model. Starting with 
an identical initial condition for each run, the equili- 
brated flow behavior displays a variety of complex fluid 
behavior as the outer parameter (in this case, rotation 
rate) is changed. The wavenumber transition regimes 
(which are characterized by one or more of the follow- 
ing behaviors: wave dispersion, two-stage equilibra- 
tion, and irregular wavenumber selection) and the reg- 
ular wave regimes (which equilibrates to a state con- 
sisting of only one wave and its higher harmonics) are 
identified. Due to the “quiescent-start” initial condi- 
tion, the hysteresis region in Lu et al. (1994) is not 
well defined. Nevertheless, the presence of multiple so- 
lutions for certain parameter settings is evident. 

The purpose of this study is to study wavenumber 
selection and hysteresis in an unstable baroclinic flow 
using a high-resolution spectral model, which includes 
waves both longer and shorter than the linearly most 
unstable wave to allow for both upscale and downscale 
energy transfers. Section 2 describes the model setup 
and its linear properties. Section 3 describes the spec- 
tral numerical technique used in solving the governing 


equations. The effects of linearly stable modes on the 
finite-amplitude dynamics for a fixed parameter setting 
are discussed in section 4. The results are used as a 
guide for determining the appropriate truncation level 
for the study of wavenumber transition and hysteresis. 
The nonlinear solutions as a function of basic flow in- 
stability are presented in section 5. Finally, section 6 
concludes the description of the numerical experiments 
and provides further remarks. 

2. The model and linear analysis 

The model used here is identical to that in Chou and 
Loesch ( 1986b, hereafter referred to as CL) — that is, 
a quasigeostrophic model based on Eady (1949) with 
the addition of unequal Ekman dissipation at the top 
and bottom rigid boundaries. A continuously stratified 
fluid, rotating on an /plane, is contained in a channel, 
periodic in the zonal direction and bounded by vertical 
walls in the north and south. The basic state is a zonal 
flow of constant vertical shear — that is, \z in nondi- 
mensional form, where X is the vertical shear and z is 
the height. This basic state yields a conservation of qua- 
sigeostrophic potential vorticity in the interior, while 
the vertical velocity pumped out of the Ekman layers 
specifies the boundary conditions at the top and bottom 
of the channel. The governing equation and boundary 
conditions for the nondimensional perturbation stream- 
function <i> are 

S7 2 <p + S~'(f) zz = 0, (2.1a) 

<£, = ^, = 0, at y = 0, 1 (2.1b) 

<f> a = -kz<fi xz + k<p x + j^J SV 2 <j> - cf> z ), 

at z = (2.1c) 

where the overbar represents a zonal average, J(g, h) 
= gjiy - gyh x is the Jacobian describing the nonlinear 
advections of the flow, S is the (constant) stratification 
parameter, and r and r, are the Ekman dissipation as- 
sociated with the lower and upper boundary, respec- 
tively (refer to CL for more details concerning the 
model). Note that the nonlinearity enters the problem 
only through the boundary conditions (2.1c). 

The linear eigenvalue problem has been solved and 
discussed in detail in CL. Here, only the properties per- 
tinent to the current study are presented. The wave so- 
lution is sought in the form of 

<t> = F(z)e ik(x ~ c,) sin/J7ry, (2.2) 

where k is the zonal wavenumber, n is an integer de- 
noting the meridional mode, F(z) is the vertical struc- 
ture, and c is the (complex) phase speed. The linear 
growth rate &c, as a function of wavenumber and ver- 
tical shear for the gravest meridional mode n = 1 is 
shown in Fig. 1 for S = 0.1, r = 1, and r, = 0.1. The 
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flow becomes unstable when X exceeds the critical 
value of K-.min = 0.587, where the linearly most unsta- 
ble wave is k = 3. The sideband waves become unstable 
when the shear is slightly increased above rain ; that 
is, the shears for k = 4 and 2 are 0.630 and 0.645, 
respectively. A shortwave cutoff is present at k = 7.72 
beyond which the waves are always stable. For the sec- 
ond meridional mode n = 2, only wavenumbers 1-5 
are unstable, while the critical wave remains k = 3 with 
a critical shear of 2.18. For even higher meridional 
modes n 5* 3, all waves are stable. The critical shears 
for the unstable modes are listed in Table 1 for S = 0. 1 , 
r = 1, and r, =0.1. The confinement of instability to 
long waves and low meridional modes is due to the 
constraint imposed on the vertical e-folding distance, 
which is inversely proportional to [S(fc 2 + n 2 n 2 )]' 12 
[see CL, (2.8)]. 

Close examination of Fig. 1 reveals that as the flow 
becomes more supercritical the fastest growing wave 
gradually shifts to a shorter wave; for example, at X 
~ 1.5 the fastest growing wave becomes k = 4. As 
mentioned in the introduction, this poses certain im- 
portant questions: Which wave dominates at finite am- 
plitude as the flow becomes more supercritical? Will 
the dominant wave shift to a shorter wave as predicted 
by the linear theory? Or, will it shift to a longer wave 
as observed in the annulus experiment and weakly non- 
linear theory? 

Another question arising from the linear analysis is 
the effect of stable wave components in finite-ampli- 
tude equilibration. As the unstable waves grow, the 
self-interaction of these waves generates a correction 
to the mean state which, in turn, modifies the stability 
characteristics of the zonal flow. Also, the wave-wave 
interactions among the growing waves may have a pos- 
itive contribution to the nonlinear growth of certain 
wave scales. The combined quasi-linear and nonlinear 
effects may then destabilize the linearly stable waves. 
This usually occurs during the initial spinup when the 
instability is strong and the feedback is intense. These 
stable waves may eventually decay to zero as the sys- 
tem reaches its final equilibrium state. However, re- 
gardless of how ephemeral and subtle these waves are, 
they may play a profound role in determining the final 
wave state since, as will be shown later, the final wave 
state depends not only on the physical parameters of 
the model but also on how the solution is approached 
in parameter space. Also, as indicated by Cehelsky and 
Tung (1987), the presence of small-scale waves pro- 
vides a path for proper energy and vorticity cascades 
to the small scales where significant dissipative sinks 
are present. For the current parameter setting only 12 
unstable wave modes were found (k = 1-7 for n = 1 
and k = 1-5 for n = 2). It is not clear whether it is 
adequate to include only the unstable modes in the 
model, even just to qualitatively describe its finite-am- 
plitude behavior. This will be examined by varying the 


Table 1 . Linear critical shear for the unstable wave components 
for 5 = 0.1, r = 1, and r, = 0.1. 


k 

1 

n 

2 

1 

1.031 


4.694 

2 

0.645 


2.652 

3 

0.587 


2.182 

4 

0.630 


2.287 

5 

0.747 


3.666 

6 

0.976 


stable 

7 

1.596 


stable 


zonal and meridional truncation levels, while keeping 
the model parameters at fixed values. 

3. The spectral numerical model 

As in CL, a spectral numerical method based on Fou- 
rier expansion is used to solve (2.1 ). The perturbation 
streamfunction, which satisfies the boundary conditions 
(2.1b), is represented by two truncated Fourier series 

M N 

<t> = X L CmU, t)e imkx Sin/J7ry 

m=l 1 

L 

+ c.c. + Y, C,(z, t) cosl-rcy, (3.1 ) 
/= i 

where the first series represents a wavelike stream- 
function, the second series represents the correction to 
the mean flow, and c.c. represents the complex conju- 
gate of the previous term. The vertical structures of the 
complex wave component C m „ and the zonally mean 
component C, are solved directly from (2.1a) to yield 

mn A rtm (r) cosh2^,„,- -F b nl/l ( f ) sinf^/u^z 
C, = M,(t) sinh2a,z + N,(t) cosh2a/Z, 

where 

^ = 0.5 [S(m 2 k 2 + n 2 7r 2 )]' /2 and a, = 0.5S U2 Itt. 

The resulting 2(M X N + L) coupled first-order or- 
dinary differential equations are integrated numerically 
in time using a fourth-order Runge-Kutta scheme. A 
transform method, which is based on the fast Fourier 
transform to compute the Jacobian terms in (2.1c), has 
been developed to replace the interaction coefficient 
method used in previous studies. This new algorithm 
increases the efficiency of the calculation and enables 
us to include enough wave harmonics in both the zonal 
and meridional directions to critically examine the ef- 
fect of stable short waves in wave equilibration. 

Unlike in the majority of the previous baroclinic in- 
stability studies by the author (Chou and Loesch 
1986a,b; 1991 ), where the primary interest was on the 
linearly most unstable wave (Jfc = 3), in the present 
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study the fundamental wavenumber is chosen as k = 1 
to include the long waves and allow for both upscale 
and downscale energy transfers from the most unstable 
wave. The basic parameters chosen for the present 
study are fixed at S = 0.1, r = 1, and r, = 0.1, which 
are close to the atmospheric values according to CL. 
The nondimensional time step used is At = 0.02, which 
has been proven to be adequate for the present study. 
The remaining unspecified variable in (2.1) is the ba- 
sic-state shear X, which measures the baroclinicity of 
the system. Here, we define the supercriticality as the 
amount of shear exceeding the minimum critical shear; 
that is, 

A X X c min • 

Due to the meridional structure of the perturbation 
imposed by the sidewall boundary conditions, the sys- 
tem can be divided into m + n = even and m + n 
= odd subsets, which interact with each other indirectly 
through the mean flow. If the perturbation contains only 
one of the subsets initially, the other subset will remain 
zero for all time. To ensure that each subset has an 
equal opportunity to grow, the initial conditions are 
chosen such that there is at least one nonzero wave 
amplitude in each subset at the beginning of each in- 
tegration. 

4. Effect of stable wave components 

Intuition may suggest that the inclusion of all unsta- 
ble modes in the wave spectrum is sufficient to quali- 
tatively describe the nonlinear wave evolution and, as 
more modes are included, the accuracy of the approx- 
imation may improve. To test this conjecture, a series 
of integrations is conducted using various truncation 
levels. The supercriticality is fixed at A = 1.2 (i.e., X 
= 1.787), which is unstable to the wave components 
with fundamental meridional mode n = 1, but is stable 
to the wave components with higher meridional modes 
(see Table 1). The initial conditions are A 2] (0) 
= A 31 (0) = 0.001. 

When the wave field consists of only the (linearly) 
unstable components and the mean-flow correction is 
limited to the gravest meridional mode, that is, (M , N, 
L) - (7, 1, 1 ), the solution grows out of bound due to 
the insufficient mean-flow representation to stabilize 
the system. Including more meridional harmonics in 
the mean field results in a damped vacillation consisting 
of a dominant wavenumber 3 and a weaker wavenum- 
ber 2. This is shown in Fig. 2a for the (7, 1, 3) trun- 
cation. Note that waves that do not exist initially are 
not generated since the restriction to a single meridional 
mode in the wave field precludes wave -wave inter- 
actions. The nonlinear interactions in this case are lim- 
ited to the self-interactions of wavenumbers 2 and 3, 
which contribute to the correction of the mean flow, 
and the interactions of the mean flow and wavenumbers 
2 and 3, which contribute to the equilibration of each 


wave. Because wavenumber 3 has a larger growth rate, 
it extracts most of the available potential energy from 
the zonal flow and contributes to the bulk of the mean- 
flow correction. In fact, the time evolution of wave- 
number 3 is almost identical to the case in which wave- 
number 3 is the sole wave in the system. The solution 
essentially remains the same for L 3, indicating a 
fast-converging series of the mean flow when the wave 
field is limited to a single meridional mode. 

The effects of stable wave components on wave 
equilibration are examined by systematically increas- 
ing the truncation level in the wave field, first in the 
meridional direction and then in the zonal direction. 
When higher meridional harmonics are included in the 
wave field, the wave-wave interaction becomes pos- 
sible and energy can be transferred in both upscale and 
downscale directions. The energy transfer is most vig- 
orous during the initial adjustment period when the in- 
stability is strongest. As shown in Fig. 2b for the (7, 2, 
3) truncation, all wave modes are generated through 
wave-wave interactions when the second meridional 
harmonics are included in the wave field. Here, wave- 
number 3 reaches a steady state through a damped vac- 
illation similar to the (7, 1,3) truncation. However, 
wavenumber 2, together with all other waves, grows 
initially but eventually decays to zero. The result re- 
mains the same when more mean field harmonics are 
included. As shown in Fig. 2c, the solution for the (7, 
3, 3) truncation changes drastically from the lower- 
order solutions. The interaction is strong at the initial 
stage, and the dominant wavenumber switches to wave- 
number 2 after about 100 time units. The system even- 
tually reaches a multiwave state, which vacillates cha- 
otically. When the wave field is better resolved, most 
of the chaos disappears as is shown in Fig. 2d for the 
(7, 5, 5 ) truncation. 1 The presence of ‘ ‘spurious chaos” 
indicates that the lower-order model does not have the 
proper channels through which energy can cascade 
upscale and enstrophy can cascade downscale (Cehel- 
sky and Tung 1987). In the higher-resolution model, 
the chaotic behavior is limited to the initial adjustment 
stage when the waves are competing for the available 
potential energy released by the instability. Afterward, 
wavenumber 2 emerges as the sole survivor of the sys- 
tem and reaches a steady state. 

A further increase in meridional resolution produces 
a similar solution with better approximation, except for 
the (7, 11, 11) truncation, 2 where a qualitatively dif- 
ferent behavior emerges. Instead of diminishing to 
zero, wavenumber 3 remains dominant for all time and 


1 Since the FFT algorithm uses only one meridional length for both 
wave and mean fields, it will be to our advantage to set N = L in 
further calculations. 

2 The seemingly bizarre number used here is due to the FFT limit 
on the total length of transformation, including the mean field, to 
2 p y5 r , where p » 1, and q, r 3= 0. 
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Fig. 2. Solutions for various meridional truncations when the zonal truncation is kept at M = 7: (a) N = 1, L = 3; (b) N = 2, L = 3; 
(c) A' = L = 3; (d) A' = L = 5; (e) /V = L = 11; and (f) N = L — 23. The parameter setting is S = 0. 1, r = 1, r, =0.1, and A = 1.2. 


vacillates perpetually, while all other waves reach small 
but finite amplitude. This vacillatory mixed-wave state 
is shown in Fig. 2e. It was puzzling at first that the 
solution converged to a totally different attractor as the 
truncation level was varied. However, it will become 
clear later that the unexpected change in solution char- 
acteristics was due to the fact that in certain parameter 
ranges a dynamic system can possess more than one 
stable solution, which is sensitive not only to physical 
parameters and initial conditions, but also to numerical 
techniques such as truncation level and integration 
scheme. 

In a chaotic dynamic system, a slight difference 
in initial condition can lead to a profound conse- 
quence in the final solution (Lorenz 1962). The 


sensitivity to initial condition can be applied to 
the current study as follows. When started with 
an identical initial condition, the difference in 
truncation level creates a different solution imme- 
diately following the first time step. Its effect is 
continuously felt at subsequent time integrations, 
contributing to the further deviation of solution tra- 
jectory in phase space. In a system where several 
stable attractors exist, the trajectory may close in 
upon one of the possible attractors and render a dif- 
ferent flow behavior when truncation level is varied. 
For reference, Fig. 2f provides the solution with me- 
ridional truncation N = L = 23, which has been 
proven to be adequate for the current zonal trunca- 
tion M = 7. 
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To investigate the effect of stable zonal wave on 
wave equilibration, the meridional truncation is fixed 
at /V = L = 23. When more zonal waves are included 
in the wave spectrum, the solution changes character 
again. This is shown in Fig. 3a for the (9, 23, 23) trun- 
cation. As in the low zonal wave case, after a period 
of adjustment wavenumber 2 emerges as the dominant 
wave. However, unlike in Fig. 2f, it vacillates forever, 
while all other waves gradually attain finite (although 
small) amplitude, resulting in a mixed-wave vacillatory 
state. When zonal resolution of the model is increased 
to M = 11, wavenumber 3 replaces wavenumber 2 as 
the dominant wave and vacillates perpetually. This is 
shown in Fig. 3b for the (11, 23, 23) truncation. In 
contrast to the dominant wavenumber 3 state in Fig. 2e, 
here all the other waves diminish to zero and the system 
reaches a single wave vacillatory state. A further in- 
crease in zonal resolution results in a mixed-wave vac- 
illatory state with dominant wavenumber 2, which is 
not dissimilar to Fig. 3a. The final converged solution 
for A = 1.2 is shown in Fig. 3c for the (23, 23, 23) 
truncation. 

Here, we demonstrate that in wave equilibration the 
higher stable modes (both meridional and latitudinal) 
serve not only to better approximate the model solu- 
tions, but also to determine which of the possible so- 
lutions is ultimately reached. Due to the chaotic nature 
of the transition period, a slight difference in solution 
trajectory at initial stage may lead to a completely dif- 
ferent final state in a multiple equilibria system. 

To determine the truncation level for further inves- 
tigation, several integrations with various truncation 
levels were conducted. It is found from previous cal- 
culations that for the solution to converge properly, the 
meridional and zonal harmonics need to increase con- 
currently; that is, M = N = L is required. Since it is 
nearly impossible to duplicate exactly the whole wave 
history from beginning to end due to the chaotic nature 
of the transition period, a strategy was adopted that 
compares only the final equilibrated state and circum- 
vents the capricious behavior of the transition period. 
A control case was run with a high resolution, (M, N, 
JL) = (31, 31, 31), until the final equilibrium state was 
reached. The final state was then used as the initial 
condition for lower-resolution integrations. By com- 
paring the equilibrium state of the control solution and 
the lower-order solutions, the optimal truncation level 
can be determined. The lowest truncation required to 
qualitatively approximate the converged solution was 
found to be (9, 9, 9). Quantitatively, its solution is 
surprisingly accurate and lies within 2% of the control 
case. For higher-order calculations, the accuracy is in- 
creased to better than 1 % for the ( 1 7, 17, 17) truncation 
and about 0.1% for the (23, 23, 23) truncation. To en- 
sure proper convergence for higher supercritical cases, 
the (23, 23, 23) truncation was chosen for the hyster- 
esis study in the next section. This truncation has been 
proven to be more than adequate, since in the hysteresis 



Fig. 3. Solutions for various zonal truncations when the meridional 
truncations are kept at N = L = 23: (a) M = 9, (b) M = 11, and (c) 
M = 23. 


study the initial condition is never far from the con- 
verged solution due to the small increment in super- 
criticality. This truncation has a better resolution than 
that in CL, where the highest retained zonal wavenum- 
ber is 18 (for k = 3 and M = 6) with 12 meridional 
modes. 

5. Hysteresis study 

In this section, the model response is examined using 
a hysteresis approach; that is, the numerical integra- 
tions are carried out by stepwise varying one of the 
external parameters with a small increment, while using 
the equilibrated state of the previous parameter value 
as the initial state of the new run. If hysteresis exists, 
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Fig. 4. Schematic diagram outlining the stable solutions of (2.1). 
Each solution branch is labeled by an alphanumeric symbol: the first 
letter denotes a single (S) or mixed (M) wave state; the second de- 
notes a steady state (S) or vacillation (V), followed by the dominant 
wavenumber. An s or w is added on the MV-2 branches to indicate 
a strong or weak vacillation. The asterisk denotes an intermediate 
mixed-wave solution, and the arrow indicates the direction of incre- 
ment. 


the model response at a given parameter value will con- 
verge to a different fluid behavior when it is traversed 
from the opposite direction. The parameter chosen for 
this study is supercriticality A, and the numerical ex- 
periment starts from A = 0. 1 to A = 2.4 then back to 
A = 0.05. In general, the increment in A is 0.1, but is 
reduced to 0.05 when a change in solution character- 
istics is encountered to ensure that the new solution is 
genuine and not a numerical artifact caused by the im- 
balance of the flow at the initial stage. 

As mentioned earlier, the wave spectrum consists of 
two nearly independent subsets, and if either one of 
them contains no initial energy or external forcing, that 
particular subset will remain zero forever. Therefore, 
when the final state of previous A contains only one 
subset, an additional amplitude of 0.001 is added in 
either A 2 1 or A 3I at the beginning of the new integration 
to ensure the inclusion of all possible growing modes. 
It was found in the subsequent experiments that the 
additional energy spreads to the rest of the subset in a 
matter of few time steps, if the instability so warrants. 

The results of a hysteresis experiment are summa- 
rized in Fig. 4. The horizontal axis is the supercriticality 
A and the vertical axis is the final dominant wavenum- 
ber. Each solution branch is given by an alphanumeric 
symbol: the first letter denotes a single (S) or mixed 
(M) wave state, the second a steady state (S) or vac- 
illation (V), followed by a number corresponding to 
the final dominant wavenumber. An s or w is added on 
the MV-2 branches to distinguish a strong vacillation 
from a weak one. 

Figure 4 shows that at finite amplitude, instead of 
shifting to the faster growing shorter wave as predicted 


by the linear stability theory, there is a succession of 
transition to longer waves as the flow becomes more 
supercritical. An example of wavenumber transition is 
given in Fig. 5 for A = 1 .35 where the flow behavior 
changes from a single wavenumber 3 state to a multi- 
wave state with dominant wavenumber 2. Note that the 
time-averaged amplitude reached by wavenumber 2 af- 
ter transition is larger than that reached by wavenumber 
3 before transition, indicating a higher efficiency in re- 
ducing instability by the longer wave. Before the 
threshold A for wavenumber transition is reached, at 
A = 1.3 the flow goes through an intermediate state in 
which the equilibrated wave state consists of a full 
spectrum in a vacillatory manner (marked by an aster- 
isk in Fig. 4). The presence of a precursory mixed- 
wave state in the transition from a single-wave state to 
another state of a longer dominant wave has also been 
observed by Hart ( 1981 ) . However, in his study, tran- 
sition through mixed-wave state and hysteresis is mu- 
tually exclusive; that is, the wavenumber transition ei- 
ther takes place through a mixed-wave state with no 
hysteresis or occurs abruptly in a hysteresis fashion 
without going through a mixed-wave state. In the pres- 
ent study the intermediate state is part of the hysteresis 
loop. 

The reverse transition, that is, the transition to a 
shorter wave as supercriticality is decreased, takes 
place at a smaller A than the forward transition does. 




Fig. 5. Amplitude and mean shear evolutions for A = 1 .35, showing 
a dominant wavenumber 3 to 2 transition. 
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Fig. 6. Three stable solutions with dominant wavenumber 2 for A 
= 1.8: (a) the MS-2 state, (b) the MV-2w state, and (c) the MV-2s 
state. 


demonstrating a hysteresis behavior. Unique solutions 
exist only when the flow is weakly supercritical (A 
0.1) or strongly supercritical (A 3= 2.25), while 
multiple stable solutions prevail in the intermediate su- 
percritical cases. Most notably, for 0.8 « A =£ 1.3 there 
exists three stable solutions, each with dominant wave- 
number 1, 2, and 3, respectively, while for 1.7 A 
=s 1.95 there are four stable solutions: a dominant 
wavenumber 1 state and three dominant wavenumber 
2 states, each with different wave characteristics. It is 
not clear whether all of the possible solutions in the 
parameter range have been found because the MV-2w 
state would have been totally missed if the finer A in- 
crement of 0.05 were used exclusively throughout the 
study. 


The three dominant wavenumber 2 solutions are il- 
lustrated in Fig. 6 for A = 1.8, each with three leading 
wave components. All three solutions reach a multi- 
wave state, but of different wave composition: the 
steady-state solution (the MS-2 solution) contains a 
full wave spectrum, the weak vacillation solution (the 
MV-2w solution) contains only even number zonal 
waves and their meridional harmonics, and the strong 
vacillation solution (the MV-2s solution) contains even 
less wave components in which the leading four wave 
components are ( m , n) = (2, 1), (6, 1), (4, 2), and 
(2, 3). The vacillations found in the vacillatory wave- 
number 2 states are quite regular in nature. The steady- 
state solution becomes vacillatory as the threshold A 
for wavenumber transition is approached (marked by 
an asterisk in Fig. 4 at A = 1.95), whose amplitude 
variation lies between those of the strong and weak 
vacillation cases. The transition from each state to the 
dominant wavenumber 1 state occurs at a different A: 
A = 2 for the MS-2 and SV-2 solutions, and 2.25 for 
the MV-2 solution. 

Although only one dominant wavenumber 1 solution 
has been found, the MV-1 branch contains a variety of 
flow behaviors. The solution type and the leading equil- 
ibrated amplitudes with their vacillation range are il- 
lustrated in Fig. 7 for 0.8 A =£ 1 .6. In the first chaotic 

region (0.8 s; A =s 0.85), the wave state displays a 
low-frequency variability, where the long waves 
(waves 1 and 2) vacillate with longer periods than the 
shorter waves (Fig. 8). The vacillation becomes less 
vigorous as A is increased, and at A = 0.9 the vacil- 
lation of wavenumber 1 becomes so weak that it ap- 
pears to reach a steady state, while wavenumber 3 vac- 
illates regularly (the range of amplitude vacillation is 
shown in Fig. 7). (Although the terms regular vacil- 
lation and double-period vacillation are used here, they 
nevertheless contain a low-level noise.) When A is fur- 
ther increased to A = 0.95, a different flow behavior 
starts to emerge (the A2 solution in Fig. 7). Here, the 
flow behavior changes aperiodically between two dis- 



Fig. 7. Solution type and leading final amplitudes for the 
MV-1 solutions. The I-bar indicates the range of vacillation. 
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Fig. 8. Amplitude and mean shear evolutions for A = 0.85, showing 
a chaotic vacillation with two predominant frequencies. 


tinct states: a relatively regular vacillation and a chaotic 
vacillation. As illustrated in Fig. 9, the chaotic vacil- 
lation is characterized by a slightly elevated dominant 
wave amplitude, a doubled period for the shorter 
waves, and a weaker zonal shear at the center of chan- 
nel, when compared to the regular vacillation. The be- 
havior of switching aperiodically between two unstable 
attractors has also been observed by Lorenz ( 1990) in 
the form of interannual variability, and by Reinhold and 
Pierrehumbert ( 1982) in the form of multiple weather 
regimes, except here it is obtained strictly through the 
nonlinear baroclinic dynamics without the introduction 
of externally imposed forcing, such as seasonally vary- 
ing heating or topographic forcing. The switch between 
two states becomes less pronounced as A is further 
increased, and at A = 1.1 the chaotic vacillation ap- 
pears again. At A = 1.2, the chaotic vacillation coa- 
lesces to a double-period vacillation, which lasts until 
A = 1.4 before transforming to a single-period vacil- 
lation. At this point, wavenumber 3 overtakes wave- 
number 2 as the second most dominant wave, while 
wavenumber 1 still remains most dominant. At A 
= 2.35, this vacillation starts to modulate with a long- 
period (~100 time units) oscillation and the solution 
takes a form of chaotic vacillation, which becomes 
more vigorous as A is further increased. 


The MV-1 solutions obtained so far seem to indicate 
three possible routes to chaos: through an abrupt 
change in flow characteristics (from A = 0.9 to 0.85 
and A = 2.3 to 2.35); through a multiple-regime stage 
where chaotic and regular vacillations appears alter- 
natively (from A = 0.9 to 1.05); and through a period 
doubling (from A = 1.5 to 1.05). A detailed investi- 
gation of the routes to chaos requires a smaller incre- 
ment in A, which is beyond the scope of the present 
study. 

Figure 10 illustrates the time-averaged zonally mean 
shear at the center of the channel for the solution 
branches shown in Fig. 4. In general, the solution with 
a longer dominant wave equilibrates to a larger mean 
shear and prevails at larger values of A. For a given 
dominant wavenumber, the mean shear decreases as A 
is increased until a minimum is reached and then in- 
creases with A. For example, on the SS-3/SV-3 solu- 
tion branch the zonally mean shear decreases from near 
the critical shear (0.587) for A = 0.05 to a minimum 
of 0.424 at A = 0.8 and then increases to 0.507 at A 
= 1.3 before it switches to the SV-2 state. The corre- 
sponding time-averaged amplitude for the SS-3/SV-3 
branch is given in Fig. 1 1 . The time-averaged wave- 
number 3 amplitude increases rapidly with A when A 
is small, levels off when A becomes larger, and starts 
to decrease after A = 1.2. This indicates that wave- 
number 3 becomes less effective in removing the im- 
posed instability as A is increased toward the threshold 
A. Because longer waves can prevail at a higher equil- 
ibrated shear, the natural selection process favors the 
transition to a longer wave state when the short wave 
becomes incapable of neutralizing the imposed insta- 
bility. 

6. Conclusions and further remarks 

A continuously stratified quasigeostrophic model has 
been used to investigate the wavenumber transition and 
hysteresis in unstable baroclinic flows. It is found that 
as the flow becomes more supercritical, the final dom- 
inant wavenumber may not be the linearly fastest grow- 
ing wave, but will shift to a longer wave. In certain 
parameter ranges, the realized finite amplitude state is 
not uniquely determined by the external parameters, 
but also depends on the direction the solution is tra- 
versed in parameter space, thus, displaying a hysteresis 
behavior. Once a flow is established at a given point in 
the parameter space, it will remain unchanged until a 
critical value of the imposed parameter is reached for 
the onset of another state. Within the hysteresis range, 
multiple stable solutions are possible and each solution 
constitutes a branch of several intertwined hysteresis 
loops. In contrast to forced wave cases (Reinhold and 
Pierrehumbert 1982; Lorenz 1990), it is shown here 
that multiple weather regime or transitive variability 
can be found strictly in a thermally driven baroclinic 
system without the externally imposed forcing. 
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Fig. 9. Amplitude and mean shear evolutions for A = 1, showing an aperiodic switch between two distinct states. 


During the course of the investigation it is also found 
that, other than the initial condition discussed in the 
previous section, the finite-amplitude state of a baro- 
clinic flow at a given parameter setting is also sensitive 
to several other factors such as truncation level and 
numerical techniques. 

As discussed in section 4, within the hysteresis range 
the finite amplitude solution may reach a different flow 
state when a different truncation level is used. If the 
prescribed truncation level is high enough, the finite- 



Fig. 10. Time-averaged equilibrium mean shear at the center 
of the channel for the solution branches shown in Fig. 4. 


amplitude state may converge to one of the stable so- 
lutions, but if the truncation level is low, the finite- 
amplitude state may approach a spurious solution due 
to insufficient wave components for energy and enstro- 
phy cascade. A preliminary experiment, conducted 
with an M = 6 and N = L = 12 truncation, shows two 
separate hysteresis loops: a wavenumber 3 to 2 transi- 
tion loop in 0. 1 =s A 0.3 and a wavenumber 2 to 1 
transition loop in 0.8 =s A 1.75. Single-wave steady 
states are always observed on the first loop, while sin- 
gle-wave steady states prevail for A 1 .2 on the dom- 



Fig. 1 1 . Time-averaged amplitude for the dominant wavenumber 
3 solution branch. The I-bar indicates the range of vacillation. 
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inant wavenumber 2 branch and mixed-wave vacilla- 
tory states prevail in the rest of the second loop. Com- 
pared to Fig. 4, the dominant wavenumber 3 in the 
lower truncation case exists in a much smaller super- 
criticality range and only one dominant wavenumber 2 
branch was ever found so that no more than two solu- 
tions are possible at a given supercriticality. The sen- 
sitivity of solution to truncation level indicates that the 
low-resolution model should be used with caution, es- 
pecially when applied to strongly nonlinear flows. 

The presence of multiple solutions also makes the 
finite-amplitude solution sensitive to the numerical 
method used, especially when starting with arbitrary 
initial amplitudes or when the parameter setting is near 
the transition point in the regime diagram. For example, 
using a moderate truncation level (M = 6, N = L 
= 12) with initial conditions A 3 |(0) = A 2 i(0) = 0.001, 
the solution for A = 1 .2 converges to a dominant wave- 
number 2 state. However, the interaction coefficient 
method yields a damped vacillation, while the trans- 
form method yields a perpetual vacillation. The sensi- 
tivity to numerical method is also observed in a nu- 
merical model of the baroclinic annulus when different 
finite difference schemes are used (Lu et al. 1994). To 
a lesser degree, perhaps, the finite-amplitude state may 
also depend on a time step used for integration, and the 
inherent architecture and round-off error of the com- 
puter used, especially when an arbitrary initial condi- 
tion is used. 

The presence of multiple solution in a quasigeo- 
strophic system raises a question of the predictability 
of a climate system, just as the presence of chaos does 
of the predictability of the weather. The sensitive de- 
pendence of the long-term equilibrium state (i.e., cli- 
mate) on initial conditions and numerical techniques 
may pose a problem in climate prediction using a gen- 
eral circulation model (GCM), since it is not clear 
whether the GCM solutions in the parameter range of 
realistic atmospheric situation are unique. The in- 
creased sophistication in GCM may act to “lock” in 
the system to a unique solution, but this argument re- 
mains illusive and further investigation is required. The 
hysteresis approach used in the present study can be 
applied to a GCM, but it may be prohibitively expen- 
sive at this moment. A more feasible alternative is to 


increase the complexity of the current model to include 
more realistic physical mechanisms, such as seasonal 
variation, land-sea contrast, and/or topographic forc- 
ing. This approach is currently under investigation. 
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